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I. INTRODUCTION 

The interpretation of the results of ultrarelativistic 
heavy ion collisions is presently one of the most chal- 
lenging problems in theoretical nuclear physics. In these 
collisions, investigated at the Relativistic Heavy Ion Col- 
lider (RHIC) at the Brookhaven National Laboratory and 
at the Large Hadron Collider (LHC) at CERN, more than 
thousand particles are observed in central collisions. Al- 
though the multiplicity and the single particle transverse 
momentum spectra at midrapidity of the different par- 
ticle species are of interest by their own right, the pur- 
pose of the experiments is to find out whether during 
the reaction the matter has made a transition towards 
a new state of matter, a Quark-Gluon Plasma (QGP). 
This information is not directly visible in the measured 
hadron spectra and therefore theoretical approaches have 
to be employed to verify whether the measured observ- 
ables are compatible with the existence of such a QGP 
or, even more wanted, whether they can even lead to the 
conclusion that such a state is necessary to explain the 
measured quantities. 

The state of the art theoretical approaches aim at a 
complete description of the heavy ion reaction, from the 
initial separation of projectile and target up to the mo- 
menta of the finally observed particles [1, 2]. Almost all 
of the presently developed models assume that quite early 
in the heavy ion reaction a QGP in local equilibrium is es- 
tablished and that hence the further evolution of the sys- 
tem can be described by the hydrodynamical equations 
before finally hadronisation takes place. Employing the 
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equation of state, calculated by lattice gauge calculations, 
taking into account initial state fluctuations and assum- 
ing that the Cooper-Frye formalism describes hadronisa- 
tion correctly, these models have been quite successful in 
describing many of the observables. 

To bridge the yet unknown physics these models neces- 
sarily rely on a number of assumptions which are highly 
debated. They assume that very early during the interac- 
tion a local equilibrium is established which is hard to un- 
derstand in view of the magnitude of perturbative Quan- 
tum Chromo Dynamics (pQCD) cross sections. How this 
can happen is therefore not clear yet. They also assume 
that no correlations among plasma constituents are built 
up during the expansion and that therefore the hadron 
production is just determined by temperature, chemical 
potential and the collective velocity of the fluid cell at 
the moment of hadronization which is considered as in- 
stantaneously. 

Therefore it may well be that these hydrodynamical 
models don't contain the correct physics. As an example 
for the ambiguity of the theoretical interpretation of ex- 
perimental results we just mention the centrality depen- 
dence of the elliptical flow, one of the key observables, 
which is equally well reproduced in three quite different 
approaches. In approaches which use viscous hydrody- 
namics [3] this centrality dependence serves to determine 
the viscosity of the QGP and therefore of the interaction 
among the constituents of the plasma, the quarks and 
the gluons. In ideal hydrodynamical approaches with 
fluctuating initial conditions [2, 4], in which only regions 
of a high energy density form a QGP, this centrality de- 
pendence is due to the impact parameter dependence of 
the relative contributions of high energy density and of 
low energy density regions. Finally, the centrality depen- 
dence is also well described in the core-corona model [4] 
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which assumes that nucleons which suffer from only one 
hard initial scattering fragment like a proton in p-p col- 
lisions whereas the rest forms a QGP whose properties 
are impact parameter independent. 

Due to these facts it is worthwhile to test the assump- 
tions on which the application of hydrodynamical models 
is based. In order to do this one needs models which do 
not assume right from the beginning that a local equilib- 
rium is established. The Parton Hadron String Dynam- 
ics (PHSD) [5-7] is such a model which allows to study 
the plasma evolution by solving a Boltzmann type equa- 
tion. There potentials between the plasma constituents 
are chosen in such a way that the equation of state from 
lattice calculations is respected. Cross sections can be 
derived from the space like part of the interaction and 
are employed for the scattering interactions among the 
plasma constituents. In this model gluons as well as 
quarks acquire a large mass when approaching the phase 
transition. Therefore the prehadrons which are created 
at the phase transition are rather heavy. Pions and other 
light hadrons are produced by the decay of these pre- 
hadrons. Another model which allows for these studies is 
a gluonic cascades realized in the Boltzmann Approach to 
Multi Parton Scattering (BAMPS) [8]. The gluon emis- 
sions and interactions during the expansion of the QGP 
move the system toward equilibrium. 

A while ago a third approach has been advanced [9] 
which is based on the Nambu-Jona-Lasinio (NJL) La- 
grangian [10, 11]. This Lagrangian is an approximation 
to the QCD Lagrangian which respects all its symme- 
tries. In the version which includes a Polyakov loop 
(PNJL) this approach also describes the equation of state 
of the lattice QCD data [12]. It has the advantage that 
all free parameters of the Lagrangian can be determined 
by static meson properties, like meson masses and decay 
constants. It contains no explicit gluons assuming that 
the in-medium mass of the gluons is large as compared 
to the transferred momentum and therefore the inter- 
action of the quarks is effectively a contact interaction. 
The quarks interact by scalar fields and by cross sec- 
tions which can be as well derived from the Lagrangian 
[13, 14]. Mesons can be produced, even in the decon- 
fined phase, but they are unstable there and may decay. 
Only when the system approaches the cross over (for low 
chemical potential the (P)NJL Lagragian shows a cross 
over and not a phase transition) the finite width of the 
meson mass disappears and stable mesons can emerge 
from the system. The light mesons are directly produced 
by qq scattering. Similar to the PHSD approach this 
Lagrangian offers therefore the opportunity to study the 
evolution of the system from the creation of a plasma up 
to the finally observed mesons. Using a N-body molecu- 
lar dynamics approach it is possible to study correlations 
and fluctuations which are built up during the expan- 
sion phase and to investigate whether observables can be 
identified which are sensitive to them. 

The cross sections calculated in this approach are quite 
small deep inside the plasma phase but, due to s channel 



resonances, they are quite large close to the cross over 
[14] where the system behaves like a liquid. Deep in the 
plasma phase the particles have only their bare mass and 
move with a velocity which is close to the velocity of 
light. In order to study the time evolution of the N-body 
system with the NJL Lagrangian we have therefore to 
develop a molecular dynamics approach for interacting 
particles which move relativistically. Such an approach 
was advanced in the original paper of the Relativistic 
Quantum Molecular Dynamics (RQMD) [15] approach 
but has never been used in practice because of concep- 
tual and of numerical problems. Some of the conceptual 
problems are related to the choice of constraints which 
one has to impose to construct such a relativistic molec- 
ular dynamics. 

The papers which contain the mathematical tools to 
develop a relativistic molecular dynamics approach are 
widely scattered. Therefore, and in order to present a 
comprehensive approach, we will start out in the next sec- 
tion with a presentation of the formalism and its deriva- 
tion. We will describe how a relativistic molecular dy- 
namics can be developed, how one can avoid the No In- 
teraction Theorem (NIT) and how the Dirac approach for 
Hamiltonian system with constraints is of importance for 
the development of a relativistic dynamics. Finally we 
present the formalism which is used. We discuss the con- 
straints and their consequences for the dynamics. The 
third section presents the NJL model as far as it is nec- 
essary to understand our approach, in particular we dis- 
cuss how masses and cross sections can be calculated. 
The fourth section is devoted to the details of the numer- 
ical realization of the approach. In the fifth section we 
present how we validated the program and some results. 
Finally, in the sixth section, we draw our conclusions. 



II. RELATIVISTIC QUANTUM MOLECULAR 
DYNAMICS 

A. Molecular Dynamics 

1. 1-body Classical Molecular Dynamics 

In the classical molecular dynamics approach particles 
are moving under the mutual influence of forces. The goal 
is to describe the trajectories of these particles in phase 
space (q^(t), Pi(t)). Knowing the phase space point for 
a given initial condition (qi(0), p^(0)) and the Hamilto- 
nian 1-L we can predict the phase space points at any 
given moment t and can calculate the value of each ob- 
servable which is defined on the classical phase space. 
The trajectory may depend in a very sensible way on the 
initial condition and may therefore become chaotic. Such 
systems are, however, not of interest here. 

We start the discussion by providing the formalism. 
We employ the Hamilton- Jacobi approach to formulate 
the motion of one particle in phase space. The equation 
of motion for an observable A, defined on the classical 
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phase space, A(q, p, t), where q, p, t are the independent 
variables, is given by 

d , OA dAdq dAdp 

The Hamilton-Jacobi equations, which present the equa- 
tions of motion of the phase space coordinates q and p 
in time can be obtained by a variational principle 



dq _ dU 
dt ~~ dp 

dp _ _<m 

dt ~ dq' 



(2) 



where %{q, p) is the Hamiltonian of the system. We can 
bring eq. (1) into the form 



dA_dA 
~dl ~ ~di 



{A,H}. 



(3) 



{^4, B} is the Poisson's bracket of A and B defined for iV 
particles as 



JV 



,, D1 ^ OA OB OA OB 

{A,B} = \—- 7^7^- 4 

t-- 1 dq k dp k op k dq k 



In the special case where A does not explicitly depend 
on time we find 



dA 
~dt 



{A,H}. 



(5) 



If we replace A by either q or p we recover the Hamilton- 
Jacobi equations, eq. (2) 



dq 

dt 
dp 

dt 



on 

Op 

<m 

dq 



(6) 



For a given initial condition (qo,Po) these equations can 
be solved, analytically or at least numerically, and we 
obtain the desired trajectory of the particle in phase 
space. For the later discussion it is important to note 
that eqs. (6) are the differential equation for the trajec- 
tory on which the energy T~L is conserved. 



2. N-body (Quantum) Molecular Dynamics 

This approach can be easily extended towards several 
mutually interacting particles and also towards Quantum 
Mechanics. The starting point for finding the time evo- 
lution of a classical N-body system is the N-body Hamil- 
tonian 



N 



N 



i i^j 



(7) 



where V(qi, Pi, q J7 Pj) is the two body potential be- 
tween particle i and j. For a given initial conditions 
(q_i(t = 0) , Pi(t = 0)) the Hamilton-Jacobi equations (2) 
can be solved analytically or numerically. This approach 
has been extended in the eighties toward the Quantum 
Molecular Dynamics (QMD) approach, a theory which 
has been successfully applied to simulate heavy ion reac- 
tions in the energy range between 50 AMeV < Ekin < 2 
AGeV [16]. This approach allowed to clarify the origin of 
multifragmentation [17], the production of mesons close 
to threshold [18] and the equation of state of hadronic 
matter [19] well above normal nuclear matter density. It 
is based on a time dependent version of the Ritz varia- 
tional principle and starts out from a trial wave function 
of a Gaussian form. The Wigner density of this trial wave 
function is defined on the phase space and has the form 



(q,-q?W) 2 



(8) 



•exp(-(p,-p^)) 2 L) 



Assuming now that the wave function of the N-body sys- 
tem is a product of the single particle wave functions and 
that the centroids q?(t),p£(t) of the Gaussians depend on 
time whereas the width is constant, the variational prin- 
ciple gives the following equations of motion 



dq? 

dt 

dp? 

~df 



d(H) 
dp? 
d(U) 
dq? ' 



(9) 



with (H) being the expectation value of the Hamiltonian 
with respect to the trial wave function. For details we 
refer to [16]. 



B. Relativist ic phase space and transformations 
between inertial systems 

1. Minkowski phase space 

One may have the idea that the equations of motion 
for relativistic particles can be obtained by replacing in 
eq. (6) the three dimensional vectors r and p by four 
dimensional vectors <2 M and This is, however, not 
true: 

a) replacing in eq. (2) q by g M and p by p M one finds 
an equation which is not covariant because T~L is the 
zero component of the energy-momentum 4- vector. 

b) eq. (2) contains a derivative with respect to the 
time t. In a relativistic theory the time is just the 
zero component of the space-time 4-vector. For TV 
particles we have furthermore N different times and 
it is not evident how these times are related to vari- 
able t of eq. (2). 
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c) these equations describe the motion of particles in 
an 8 dimensional phase space in which neither the 
energy is conserved nor the times of the different 
particles are synchronized. In a molecular dynam- 
ics approach we are interested to obtain physical 
trajectories in a 6N + 1 dimensional phase space 
( ( li( r )? Pi( r ))? i- e - world lines of the particles, and 
we want to know at which position in coordinate 
and momentum space the particle is located for 
a given value of the time evolution parameter r 
(whose nature will be discussed later). 

Thus a Hamiltonian in a nonrelativistic sense (the to- 
tal energy of the system) has no place in a relativistic 
approach. If we talk later of an "Hamiltonian" in rela- 
tivistic dynamics and of time evolution equations which 
have a form similar to eq. (2) the meaning of the differ- 
ent terms in this equation will be completely different as 
compared to that in a non-relativistic theory. 

We start out from the 4-position and 4-momentum co- 
ordinates (q%,p%) as canonical variables which obey 



{<£,<£} = {p£,p v b } = o 



(10) 



with g^ v being the Minkowski metric with the diagonal 
{1,-1,-1,-1} and zero otherwise. Here we have intro- 
duced the Poisson brackets for 4- vectors 



N 



dA dB dA dB 



k=l 



dpkfi dp% dq k ^ ' 



(11) 



Because in a dynamical system and p^ depend on the 
time evolution parameters r, these quantities have to be 
taken at equal r. 



with Auj being small. The invariance of the scalar prod- 
uct of 4-vectors under a Lorentz transformation can be 
expressed as 



(15) 



and hence 



g ap = A£s M „A£ 

= gia ,W + Au£W p +Au, v p ) (16) 
= g ap + Aoj ap + Aoj p<7 + 0(Aoj 2 ). 

Consequently Auj pu has to be antisymmetric. There are 
6 independent elements which satisfy 



Auj ap = -Auj pa . 



(17) 



In matrix form we can write the infinitesimal Lorentz 
transformation as (the factor ^ is convention in order to 
obtain the standard definition of the angular momentum 

Ji = \e ijk Mi k ) 



A(Aw„„) = 1 - X -Au pv M pv 

= 1 - ]^Au} pv {q p d u - q u d p ) 



(18) 



where M^ v = —M UfI are the generators of the Lorentz 
group and we find (compare eq. (14)) 

A(Auj^)q a = q a - ^(Ac*v - Aujp^q^ 
= q° + AuX- 



2. Poincare Group and Algebra 

Relativistic theories have to be invariant under Lorentz 
transformations A and space-time translations a. Both 
transformations form the Poincare group with the group 
element i?(A, a). It consists of all transformations of the 
form 



Similar for the infinitesimal translation 
T(Aa) = 1 - iAa^P^ 



we find: 



= 1 + Aa^d u 



q ,v = T(Aa)q u = q v + Aa v 



(20) 



(21) 



R(A, a) = A^q v + 



(12) 



which leave the scalar product between two 4-vectors un- 
changed 



Q uQ 



(13) 



with q^ = (t, q) and = g^ v q v '. 

The algebra associated with the continuous symmetry 
group is given by the algebra of the generators of in- 
finitesimal transformations. Finite transformations can 
be built with help of the infinitesimal ones. We start 
out with a Lorentz transformation which differs only in- 
finitesimally from the neutral element R(t,0) 



= 11 + Au 



(14) 



If the system is composed of several particles we find for 
the generators 



• for the translation group 

N 



(22) 



• and for the Lorentz group (SL(n = 2, C) — » dim = 
2(n 2 - 1) = 6) : 



N 



(23) 



k=l 
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These 10 generators respect the algebra of the group 
which is called Poincare algebra : 



{M MI/ , Pp} = g^Py - g up Pp 
{Mp V , M pa } = gp P M va - gp a M V p 



(24) 



gu P Mp a +g V(T Mpp. 



This can be directly verified by going back to the defi- 
nition, eqs. (18) and (20), and calculation the brackets. 
The generators Mp V and Pp do not commute. Physically, 
that comes from the fact that there is a length contrac- 
tion in the Lorentz boost (and a time dilatation). 

The generator of a Poincare transformation is given by 
the combination of that of the Lorentz transformation 
and of a translation 



G 



uj^Mp V -a»Pp. 



(25) 



If two inertial frames O and O' are connected by an in- 
finitesimal element of the Poincare group, i?(A, a), then 
the space time coordinates of the same event in O and 
O' are related by 



{q,G} 
= A^q" + a". 



(26) 



3. Reduction of the dimension of the phase space 

Relativistic theories are based on 4-vectors whose 
transformation between two inertial systems is given by 
elements of the Poincare group. This has as consequence 
that the phase space of a N particle system has not any- 
more 6N dimensions as in nonrelativistic dynamics but 
87V. World lines are given by (q*(r), P;(r)) and there- 
fore physical trajectories (position and momentum of the 
particles as a function of the time r) have 6N + 1 dimen- 
sions. Thus we need constraints to reduce the number of 
degrees of freedom in the relativistic phase space. After 
a introduction to the SN dimensional phase space and 
to the Poincare group and algebra we will illustrate the 
reduction of the degrees of freedom first for the example 
of 1 free particle and then we extend systematically the 
approach to N-interacting particles. 



C. Prom 1 to N-body relativistic system 

1. The case of 1 free particle 

We start with the most simple case of 1 free particle 
[20]. The Hamilton's equations for the time evolution 
of a nonrelativistic particle determines that trajectory in 
the phase space for which the energy is conserved. This 



suggests to define a constraint, the mass shell constraint, 
which reads for a noninteracting particle as 



K 



m 



0. 



(27) 



It reduces the phase space from 8 to 7 dimensions by 
relating the energy of the particle with its 3-momentum 
(and also with its position if we include a potential). It 
defines therefore the 7 dimensional subspace E of the 8 
dimensional phase space on which this condition is ful- 
filled. Because K is a Poincare invariant quantity, we 
find 



{K,M^} = 0, {K,Pp} = 0. 



(28) 



Of course the 7 dimensional phase space region E is also 
Poincare invariant 



R(A,a)E = E. 



(29) 



The trajectory in phase space on which this constraint is 
satisfied is given by the solution of 



Mr) 

dr 
dp(r) 

dr 



\{q(r),K} 
\{p(r),K} 



(30) 



with the initial condition q(0) = qo and p(0) = Po- ^ 
is a free parameter. In order to associate to each value 
of r one point in phase space (q(r),p(r)) or in other 
words in order to create a worldline a second constraint, 
x(^ /i ,p M ,r) = 0, has to be employed which fixes A. It 
relates the time q° of the particle with a Lorentz invariant 
system time r. This time constraint x nas been chosen 
quite differently in the literature, giving quite different 
time evolution equations. The subspace we are interested 
in is determined by a conserved x an d K constraint. This 
is expressed by 



dX dx 



dr dr 
This equation determines A 

dx 



A{ X (r),K} = 0. 



X = -f r {x,Ky\ 



(31) 



(32) 



A depends therefore on the choice of the constraint x- 
Formally we can define 



Z = \K = -^{ x ,K}- 1 K 



(33) 



and obtain a time evolution equation for a phase space 
function / 



dr dr 

= d J- 

dr 



if, Z) 



(34) 
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which is formally identical with the nonrelativistic evo- 
lution eq. (3) but Z is not the classical Hamiltonian, it 
is given by eq. (33). 

How to treat a Hamilton system with constraints has 
been developed by Dirac [21]. To determine the time 
evolution for any function of the phase space variables 
along the trajectory determined by the two constraints 
01 = K and fa = X is given by the Dirac bracket which 
is defined as 

{A B} D = {A, B} - {A, MCMiB}, (35) 

with the matrix of these constraints 

Cy 1 = ( 36 ) 

On the hypersurface, where the constraints are fulfilled, 
Dirac brackets and Poisson brackets agree. Dirac intro- 
duced the symbol « to describe {A,B}d « {A, B} that 
two functions are identical at the subspace which is de- 
fined by the constraints. For our example we find 



{A,B} D = {A,B} 



{A,K}{ X ,B} 

{K,X} 
{A, X }{K,B} 
{X,K} • 



(37) 



The Dirac brackets of the 10 generators of the Poincare 
group yield the same result as the Poisson brackets, eq. 
(24), because K commutes with them. Therefore we can 
use also the Dirac brackets to construct a transforma- 
tion between the two inertial systems O and O' . This 
transformation we call i?*(A, a). Both transformations, 
i?(A, a) as well as i?*(A, a), map therefore E to E but 
i2*(A, a) and i?(A, a) map the same point in O to differ- 
ent points in O' . Because {G,x}d = 0, x is unchanged 
under a transformation i?*(A, a) and the Dirac brack- 
ets transform a phase space point on O to a phase space 
point on O' which has the same value of r. For the trans- 
formation using Poisson brackets this is generally not the 
case. Therefore the Dirac brackets are the proper tool to 
determine world lines in the two inertial frames [22]. Us- 
ing eq. (26) and replacing the Poisson bracket {•, •} by 
the Dirac bracket {•, -}d we find the canonical transfor- 
mation between two inertial systems 



r(T) = q»(T) + {q»(T),G} D 



(38) 



where G is given by eq. (25). If we use the Poisson 
brackets for the transformation between the two inertial 
systems we obtain the geometrical transformation 

<?">) = g"(r + At) + (r + At), G} 

dq» (39) 
« g "(T) + A-Ar + { g »,G}. 

Applying the general eq. (37) we can relate {q fl ,G}D 
and {V,G} ({K,G} = 0) : 

{«", G} D = {?, G} - {e, ^ { ^ G) - (40) 



Using furthermore the time evolution eq. (34) we find 

d X {q»,K} 



dq» 

dr dr {x,K} 
and therefore eq. (40) can be rewritten in the form 



(41) 



W' - G} D = W'. G} - {x, G} ('g) ^. (42) 



Consequently, if we can ensure that 

'dx 



{X,G} 



dr 



= At 



(43) 



the transformation between two inertial systems using 
Dirac brackets (the canonical transformation, eq. (38)) 
becomes identical to that using Poisson brackets (the ge- 
ometrical transformation eq. (39)). If the condition (43) 
is fulfilled the world lines of particles remain the same un- 
der the two transformations. They are therefore frame 
independent. This requirement of the frame indepen- 
dence of the trajectories is called the World Line Condi- 
tion (WLC). 

The constraint x determines the time evolution of the 
system. If we impose the constraint x = Q° — t — [20] 
we find for the time evolution of eq. (41) 



dqv 
d7 



\{q^K) 



(44) 



whereas for the condition x 
obtain 



XfiP^ — nfiT = [23] we 



dq^ 
d7 



\{q^K} 



rn 



(45) 



In both cases we have dp^/dr = compatible with the 
fact that we have a single free particle. The different time 
evolution equations remind us that r is a parameter in- 
troduced by the constraint x an d not an independently 
defined time. Thus the time evolution of a relativistic 
system is only determined after the constraint x is i m_ 
posed. Different choices of the constraint yield a different 
time evolution of the system. 

Now we come back to the question whether we can find 
a time evolution equations for a relativistic system which 
resembles the form of the time evolution equations found 
in non-relativistic systems 



dqv 

~d7 
dp» 

"d7 



= {q",H} D = \{q*,K} 
= {p^U} D = \{p^K}. 



(46) 



where H ^ XK = Z. This can be done using the method 
described in [24]. Anyway it is always possible to find 
back such a definition using the classical limit of the rel- 
ativistic equations of motion. 
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2. Extension to 2 interacting particles 



yields 



The above discussed construction of world lines on 
which a particle moves independent of the chosen ref- 
erence system has been extended to a larger number of 
particles in [20, 23, 24]. For a system with two interact- 
ing particles [22] the mass shell constraints eq. (27) have 
the form 



K 2 V2P2V 



m 



V = 

V = 



in order to have reference frame independent world lines. 
In addition, they have to be first class constraints in the 
notation of Dirac [21] : 



{K U K 2 } = 2 



d 



d 



V = 0, 



(48) 



a condition which can be fulfilled if the potential V de- 
pends on q£ [20] which is the part of = q± — q% which 
is transverse with respect to the centre of mass motion 
p\i — + anc [ which is defined as 



P 2 



(49) 



Poincare transformations map the 77V dimensional phase 
space, on which the constraints eq. (47) are fulfilled, on 
itself. The evolution equations can be extended to 



dr 
dr 



v 1 {q?,K 1 } + v 2 {q?,K 2 } 
v 1 {p?,K 1 } + v 2 { P f *,K 2 } 



(50) 



with arbitrary parameters v\ and v^- For an interacting 
system the time components of particles q® become con- 
nected by the potential term and consequently the spatial 
position of each of the particle qf depends on both times 
q® and q®. This does not correspond to a world line but 
to a sheet. To obtain a world line we have to synchronize 
first the times of both particles by a constraint without 
any parameter 



(51) 



and finally to connect the synchronized times to a clock 
time r 



X2 (01,02, Pi ,P2,t) = 



(52) 



with the property det {Ki,Xj} 0- H {KiiXj} = we 
cannot assign to each point on the trajectory uniquely 
a value of the parameter r. With these two additional 
constraints which reduce the 77V dimensional phase space 
to a 67V dimensional phase space with a parameter r so 
effectively to a 67V + 1 dimensional phase space [25]). 
Condition 52 allows for fixing the free parameters V{ of 
eq. 50 (see eq. 34). 



dX2(r) _ dx2(r) 



dr 



dr 



+ { X2 (r),i^ = (53) 



-{X2(r),Kj 



-i gXg (r) 
dr ' 



(54) 



Consequently, the general evolution equation for a phase 
space function / is 



(47) with 



d/ _df d X 2 { . „ x 



Sij = {Xj,Ki} 1 - 



(55) 



(56) 



As in the one particle case this requires that the two 
transformations between the inertial systems, the one ex- 
pressed by Dirac brackets (canonical transformation) 



q'?(T) = q?(T) + {q?(T),G} D 



(57) 



and the one expressed by Poisson brackets (geometrical 
transformation) 



q'i(r) = «T(t + ArO + {<f{r + An), G} 
«gr(r) + ^Ar, + {<(r),G} 



(58) 



lead to points on the same world lines. Employing the 
Dirac brackets and taking advantage of {Ki,G} = we 
find 



dgf(r) 
dr 



= {qt(r),K j }S lj ^Ar i , 



(59) 



or if {q? (r), Kj} ^ and Sij + 

te,G} = ^Ar, 



(60) 



In reality the last equation poses two conditions: Ari = 
Ar2 and that xi which does not depend on r is Poincare 
invariant to fulfill = 0. These conditions cannot 

be fulfilled by every choice of Xi- Indeed, if we relate 
the r to the fourth component of g^, the instant form of 
Dirac [21], 



Xi = ^i°-<Z2°) = 

X2 = \{qU4)~r = ^ 



(61) 



we recover {xu G} =^ and hence the no go theorem 
stating that relativistic molecular dynamics can only be 
formulated for non interacting particles [20]. If, on the 
other hand, the x% are defined kinematically as 



Xi = = 



X2 = ^ + q^-T = o 



(62) 
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with = P^/VP^, which gives = (1, 0) in the centre 
of mass of two particles, the world line condition can be 
fulfilled [20] by setting An = -{x2,G}. With the latter 
time constraints eq. (62), we can compute the Vi (eq. 
(54)). We start out from the matrix of constraints : 



{Xi.tfi} {X1,K 2 } 
{X2,K l } { X 2,K 2 } 



(63) 



which can be inverted 



Si 



"(2 



(2 pj^) 
(2 P^) 



(64) 



The parameter V{ becomes 

Wl = (2K^)- 1C = 
v 2 = (2p^)- 1 = s 



2E X 
1 

2~E 2 



(65) 



and we obtain finally the equations of motion for 2 inter- 
acting particles in their center of mass 



dr 
dr 



1 dV(q T ) 
^2E k dq i}l 



(66) 



We can easily see that the classical non-relativistic limit 
of these equations gives the same result as QMD by tak- 
ing p <C m. 

In this example of two interacting particles we can also 
address another problem which is the separability of clus- 
ters. In contradiction to nonrelativistic dynamics this 
separability is not trivially fulfilled by taking a vanish- 
ing potential for large distances because the potential 
enters the constraint matrix which determines the time 
evolution. Cluster separability means that we have the 
equations of motion of two free particles if the distance 
between them is large. 



3. Extention of the formalism 

In order to define world lines for relativistic particles 
we have to reduce the relativistic phase space from 87V 
dimensions to 67V +1 dimensions. This is done by relating 
the zero component of the space time and of the energy 
momentum 4- vector, to the vector parts of these 

4- vectors with help of constraints. In addition we have to 
introduce a system time r which describes the evolution 
of the system. 

For the two body case, for which {^, Kj} = 0, the for- 
malism has been developed in the last subsection. Here 
we extend the formalism to TV > 2 where {Ki,Kj} may 



be different from 0. In this case the time evolution equa- 
tions for the 27V constraints are 



2N 



d7 = -Q^+J2 Xk ^M = 

k 



(67) 



with 



K k (qV lP v) = for 1< k < N 

Xkiq^iP^) =0 for iV+K k < 27V- 1 (68) 



Only the constraint i = 27V depend on r. Rewriting the 
last line of eq. (67) as 



2N 



(69) 



with di being a vector in which only the 2N th component 
is different from zero we find 



— —CikCLi — 



dr 



(70) 



Also this equation shows that different choices of con- 
straints will yield different values of and different 
will give a different time evolution. Therefore the rel- 
ativistic kinematics is only defined after the constraints 
are defined. This will be discussed in the next subsection. 
Defining 



2N 

5> 



k9k 



(71) 



the trajectory in phase space for which the constraints 
are fulfilled are given by 



dr 

M 
dr 



{P?(r),Z}. 



(72) 



and have therefore the form of which reminds us on the 
Hamilton-Jacobi equations. This can be easily under- 
stood by recalling that the non-relativistic equations of 
motion lead to a trajectories on which the total energy 
of the system is conserved whereas the eqs. (72) lead to 
trajectories on which the constraints 0& are conserved. 



4- N -particle system 

The N-body system is actually a trivial generalization 
of the 3-body problem. Therefore we will discuss here the 
3-body case which was studied in detail in [15, 20, 24]. 
As compared to the 2-body case we are confronted here 
with several new features: 
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• even if we know that the WLC imposes that exactly 
one of the Xi nas to depend on r, we don't know 
their definition and there are several possibilities 
that we will discuss, 

• the commutation of the on-shell mass constraints 
{Ki,Kj} which is easily fulfilled in the 2 particle 
case (eq. (48)), and which avoids that x constraints 
appear explicitly in the equations of motion, is not 
necessarily fulfilled for the 3 particle case. 

• we still have to respect the WLC and the cluster 
separability. 

We start out with the definition of some useful quan- 
tities. First of all we define the frame projectors (this 
name is used in [26]) : 



u^ 1 



(1,0,0,0) 



where p^- 



Pj and 



U 11 = 



= (1,0,0,0) 



(73) 



(74) 



where P M is the Poincare generator (22). We call the 
latter system laboratory system because in the collider 
physics the total center of mass system is identical to the 
laboratory system. 

These projectors are used to select the 4 th component 
of 4- vector in the desired frame. Here u^- is the projector 
for the centre of mass (cms) of two particles i and j, 
and is the center of mass of for the whole N-particle 
system (lab system). The derivatives of these quantities 
are : 



8ql 

dpi 



= 



u ij u ij) 



(Sik + Sjk) 



(75) 



for a more compact notation. These projectors allow for 
defining two transverse distances : 



f 1 n (T n 



and 



4r% 



i 2 



IJ J ™IJ 

2 



(79) 



=<&e 



--1% 



((qiHaUnU", 



(80) 



We notice the following properties : 

(QTij)^ = 

and 

(QTiM" = 



(81) 



(82) 



The derivatives of these transverse distances are rele- 
gated to appendix 1 a. 

We start the discussion of the 3-particles dynamics 
with the proposition of the constraints made in the Rel- 
ativistic Quantum Molecular Dynamics (RQMD) of ref. 
[15]. In this paper the authors extend the previously dis- 
cussed 2-body system to a N-body system and employ as 
mass shell constraints including the potential V 

Ki = fiViv ~ m? + VMt%) = 0. (83) 
In the same way the time constraints are defined by 



Xi 



o 



Ki<JV-l. 



(84) 



and 



We can define 



= (g 



""-«) c = 



and 



lab 



/0 
0-10 

0-1 

\o 



/0 

0-1 

0-1 

Vo o o 



(76) 



(77) 



(78) 



These constraints connect times between couples of par- 
ticles. The last constraint 



Xn 



N 







(85) 



ensures that all times are related to the time evolution 
parameter r. 

All these constraints fulfil the WLC [15]. In order to 
assure the separability of clusters two particle distances 
are weighted with the weighting function 



9ij 



L 



exp 



(86) 



There are several problems with this approach. The 
first is that the Komar-Todorov (KT) [26] condition, 
{Ki,Kj} = 0, is not fulfilled. This means that the mass 
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shell constraints of the different particles are not inde- 
pendent. 



{Ki,Kj} 



= H 



{Vi,Vj} 



(87) 



using the fact that dVi/dq^ = —dVj/dq^. We notice that 
neither V(qr) nor V(q' T ) can fulfil this condition. In the 
RQMD paper [15] it is assumed that {Ki,Kj} remains 
negligible and consequently the time constraints do not 
appear in the full equations of motion based on eq. (72) 



dr 
dr 



N 



dK k 



2N 



dXk 



k=i 

N 



dPi » k=N+l 
2N 



k=l 



d P in 
dXk 

k^r 



E A * i£r- E A * 



because if we assume {K^Kj} = 0, then \ k = for 
N + 1 < k < 2N (see eq. (70)). With this assumption 
and the time constraint eq. (84) and (85) the parameter 
A becomes 



where SNk is defined in eq. (56). The equations of mo- 
tion of [15] are then given by 



M 

dr 



N 



(90) 



>Nk~ 



k=l 



Even if we deal with three free particles and therefore 
the KT condition is trivially fulfilled, the approach of 
ref. [15] poses problems. Taking the constraints 



(91) 





= Pi 2 


— mi 2 


= 




= P2 2 


- m 2 2 


= 


K 3 


= P3 2 


- m 3 2 


= 0, 



for free particles and 

- (<7l2^12/i 



Xi 

X2 
X3 



= (^21^21/i 



<?l>13/i)/3 = 
^23^23/i)/3 = 



(92) 



(qi + Q2 + • 



Nk, 



(89) gives the matrix of constraints 



'4/3 Pi(^12 +^13)/i 

-2/3 p%ui2p 
2/3 rfUp 



-2/3 p%u 12fl 

4/3 p%(u 2 i +^23)/i 

2/3 p%U» 



-2/3 p%u 13 ^ 

-2/3 P^23/x 

2/3 p££/ M 



(93) 



whose inverse is highly non-trivial. The numerical calcu- 
lation of Xk gives non-physical trajectories with veloci- 
ties above the speed of light. Moreover, if we include the 
weight function, eq. (86), for the separability of clusters 
we encounter another numerical problem : the matrix 
inversion fails due to the fact that the numerical values 
of the matrix elements cover very many orders of magni- 
tude. This has been discussed in [27]. 

Last but not least the energy is not conserved locally in 
time (it is conserved on average over a long time) because 
of the use of qr as variable of the potential. This choice 
imposes to compute the forces in each two body centre of 
masses system. Therefore an inverse Lorentz boost has 
to be applied to transform all forces into the same com- 
mon frame. This transformation is problematic because 
the Lorentz transformation is not valid when we study 
accelerated particles and indeed creates fluctuations of 
the energy of the system. 

To avoid these problems we made a different choice of 
constraints. Instead of formulating the constraint in the 
two body systems we define them in the common center 
of mass system. These means that we have to replace the 



u^- by in the constraint formulas which then read as 

Ki =PiPi,-m 2 i +V i (q' T 2 ) =0 

Xi = -^—U v = (94) 

Xn = ^U„-t = 0. 

We cannot fulfil the KT condition using q' T (see appendix 
1 b), but the weaker condition {i^, Yli^j Kj} = 0- This 
means that the mass shell constraint of the particle i 
commutes with the sum of that of all other particles what 
is not the case in the approach of ref. [15]. Similar as in 
[15] we assume that {Ki,Kj} is negligible. 

Our choice of constraints avoids all the other problems 
of the approach of [15]. We can check that using the 
previous example of 3 free particles. Changing the time 
constraints to 

Xl = (912 + qisTU^ = 

X2 = (q 21 + q 23 rU^/3 = 1 ' 
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which gives we find the solution 

/ 4/3 -2/3 p££7 M -2/3 p£tf M 

S^ 1 = -2/3 tfU^ 4/3 -2/3 p£tf M I , (96) 

V 2/3 2/3 p^ 2/3 p^ 



(2 p^)" 1 (2 p^y 

Sij =| (2 p^)" 1 (2 p^)" 1 | . (97) 

-(2^)" 1 -(2^/7 A1 )- 1 (2^)-! 



The last column is the A parameter which has an analyt- 
ical and trivial solution 



A fe = (2 p^y 



-l lab 



1 

2Eh 



(98) 



which is in perfect agreement with the solution we found 
for the 2-particle case. The equations of motion in the 



global frame (lab), where J^f Pi = 0> are then 



dr Ei 



dpi 
dr 



N 



1 dV k (q' T ) 



(99) 



These equations conserve energy and ensure physical tra- 
jectories with velocities below the speed of light. More- 
over, the analytical solution for the is useful to avoid 
the numerical inversion of the matrix of constraints at 
each time step of the evolution, what is not possible with 
presently available computers. The equations of motion 
of eqs. (99) are finally identical to the relativistic equa- 
tions which are currently used in other approaches [28]. 

Our approach also avoids the problem of the cluster 
separability. This can easily be seen by dividing the sys- 
tem into two subsystems a and b with 



P b- 



(100) 



If we calculate the time evolution equations for the par- 
tons in each subsystem separately we obtain the same 
result as if we calculate them for the full system. This 
means that one cluster does not influence the motion of 
the other. 



III. NAMBU-JONA-LASINIO MODEL 

In this paper we study the expansion of a q/q plasma 
employing the Nambu-Jona-Lasinio (NJL) model. The 
NJL model is the simplest low energy approximation of 
QCD. It describes the interaction between two quark cur- 
rents as a point-like exchange of a perturbative gluon [29] . 
Assuming that the mass of the gluon is large as compared 
to its momentum the interaction reduces to an effective 



four point interaction and is given by 

^c-l 3 

= « ]r y.^i^ c u^s) (101) 

c=l i,j ^ ' 

where we have explicitly shown the color/Dirac a, /?, 7, S 
and flavor i,j indices. We normalize Yli=o Kxp^pa = 
Applying a Fierz transformation in color space to this 
interaction the Lagrangian separates into two pieces [30] 
: an attractive color singlet interaction between a quark 
and an antiquark (£( qq -)) and a repulsive color anti-triplet 
interaction between two quarks (£( qq ))) which disappears 
in the large N c limit. Usually a six point interaction in 
the form of the 't Hooft determinant is added (£a)) to 
break the unwanted U a (1) symmetry of the Lagrangian. 
For this study we are only interested in the color singlet 
channel (the color octet channel gives diquarks and can 
be used to study baryons [31, 32]): 



£ = £ + £ 



(102) 



£0 is the Lagrangian for a particle without interaction. 
Concentrating on the dominant scalar and pseudo scalar 
part in Dirac space we find the following explicit form of 
the Lagrangian: 



£= £ [HfW-rn f )q f 



f={u,d,s} 



(103) 



a=() 



G D [det[g/(l - j b )q f ] + det[g/(l + J^Qf}] • 



The first term is the free kinetic part, including the flavor 
dependent current quark masses which break explic- 
itly the chiral symmetry of the Lagrangian. The second 
part is the scalar /pseudoscalar interaction in the mesonic 
channel, invariant under SUa(3) ®Ua(1)- It is diagonal 
in color as the third part, the 't Hooft determinant. The 
det runs over the flavor degrees of freedom. Consequently 
the flavors become connected. Gs is the qq coupling con- 
stant and Gd the coupling constant of the 't Hooft term. 
The quarks in the NJL Lagrangian have four point in- 
teractions (with a coupling constant Gs) and 6 point 
interactions (with a coupling constant Gd)- 
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The thermodynamic properties if this Lagragian are 
summarized in ref. [11]. The NJL Lagrangian has been 
discussed im many review articles [29, 33, 34] where all 
details of this model can be found. We concentrate here 
on those quantities which enter directly in our calcula- 
tion, the quark and meson masses as well as the cross 
sections. 



A. Quark Masses 

In the NJL model the mass of a free quark of flavor i 
is given by [29] 



m % =m* i +Gs(±N c ){iTr Si) 

- G D (2N 2 C + 37V C + l)(iTr Sj)(iTr S k ), 



(104) 



where m° is the bare mass of quark, iV c is the number of 
colors, with i ^ j ^ k and 



Tr S k = Tr S k (x = 0) 



-I 



Tr S k (p) 



Tr 



1 



(27T)4 
A d A P 

(2tt) 4 " p'-m k +ie 



(105) 



r 1 d A p m k 
J (2nfE p 



with E p = \/p 2 + m\. We use here a 3-momentum cut- 
off to regularize the integrals. If the quark is brought 
into matter with a finite baryon density \i and a finite 
temperature T thermal field theory has to be employed 
and we have to replace 



P Oo, P) Pn = fan ± P) 

J (2tt)4 ZfJ (2^)3 



(106) 



where uo n = (2n + 1)ttT with n = 1,2,... are the Mat- 
subara frequencies for fermions. Hence we find for the 
propagator 



S(p) -> S(uj n ,p) 



]/> n -m + 7°/i 
j) + m 1 

2E p iuj n - (E p - fi) 
$ — m 1 



(107) 



with 



and 



2E p iu n + (Ep + /i) 
i> = j°E p - 7 p (108) 
/ = 7 °E P + 7P. (109) 




100 200 300 400 

T (MeV) 



500 



FIG. 1: Masses of ^ and s quarks as a function of T. 



This yields [29] 



Tr S k 



2 f A d 3 p m 



% J (2tt) 3 E p 



(110) 



[l-f(E p -ri-f(E p + ri] 

with the Fermi-Dirac distribution 

f(E p ± /i) = [1 + exp((^ p ± M)/T)] _1 • (HI) 

Eq. (104) allows then to calculate the quark masses 
which are displayed on Fig. 1. 

B. Meson Masses and Coupling Constants 

How appear mesons in a theory whose Lagrangian has 
only quarks as degrees of freedom ? This is shown in Fig. 
2. It displays the scattering of a quark and an antiquark 
in our theory with 4-point interactions. The left hand 
side displays the series of the exchange terms which ap- 
pear in the the random phase approximation. This series 
can be summed up. The sum is formally displayed on 
the right hand side of this figure. This sum corresponds 
in leading order of N c to the propagator of a meson with 
the proper quantum numbers. 

The central building block for the random phase ap- 
proximation is the quark-antiquark polarization propa- 
gator (Fig. 3) 



75 



fj ,J y ' 

{Tj ff ,St (p+ 7 5 (T J ) r/ S / ' (p-\q 



(112) 



where / and /' are the explicit flavour indices and tr 
refers therefore to the spinor trace only. Ti and Tj select 
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FIG. 2: Effective interaction between two quarks in Random Phase Approximation (RPA). 



-in p {iv m , k) = ^75 




*75< 



FIG. 3: The quark- ant iquark polarisation propagator for 
pseudoscalar coupling. 







FIG. 4: The meson propagator corresponds to the RPA sum. 

the appropriate flavour channel 

A 3 for tt° 



^(Ai ±i\ 2 ) for 7r + ,7r 
j=(\ 6 ±i\ 7 ) for 
^(A 4 ±iA 5 ) for K+,K~ 



(113) 



For the more complicated r\ and rf, where U p / S is not 
diagonal, we refer to [29] where this is treated in detail. 
After the traces and sums of the polarization propagator 
eq. (112) are carried out one arrives at 



\n p / s (q 2 ,m 1 ,m 2 ) = 4N c h( mi ) 



with 



h(m) 



J (27T 



4iV c /i(m 2 ) 
4:N c q 2 I 2 (q 2 ,m 1 ,m2) 

d A p 1 



(114) 



)4 p 2 



I 2 (q 2 ,m 1 ,m 2 ) 



d A p 



1 



(2tt) 4 p 1 - m 2 
1 



(115) 



(p + q) 2 



Having the polarisation propagator we can sum up the 
terms of Fig. 2. The interactions among a quark and an 



antiquark with pseudoscalar coupling in random phase 
approximation can be written as 



2iG s + 2iG s -Il P/S '2iGs 

1 

2iG s 



(116) 



1 - 2G s U p / s 



If a pseudo scalar meson with a mass M is exchanged 
between the quarks, Fig. 4, we find for the interaction 



. m -i9l g g(M) . 

175 k k 2 -M* 275 



(117) 



Eqs. (116) and (117) have the same structure and there- 
fore we can identify the exchange of a pseudoscalar meson 
with the RPA summation of qq exchanges 



2iG s 



^Trqq 



1 - 2G s n p / s k 2 -M 2 ' 



(118) 



The mass of the meson can be obtained by solving the 
equation 



2G s Il p / s \ 



k 2 =M : 



= 



(119) 



while the coupling constant g^qq can be related to the 
residue of the pole. Expanding eq. (116) around its pole 
k 2 = M 2 we find 



P/S 



2iG< 



( an 

I dk 2 



\k 2 =M 2 



1 - 2G s U p / s k 2 - M 2 

and therefore we can identify 

-1 

\k 2 =M 2 



f du p / s 

V dk 2 



2 

9nqq' 



(120) 



(121) 



For finite temperature and finite chemical potential we 
have to replace in eq. (112) the propagators S by imagi- 
nary time propagators S. 



75 (Tj ) /' jS f ' (uji + v n , p + q) . 



(122) 



The boson frequencies v n are even, v n = ±2n tt T, n = 
0,1, 2, 3,..., while the fermion frequencies uji can take 
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odd values only uo m = ±(2m + 1) tt T, m = 0, 1, 2, 3, 

So in order to find the pole mass of the pseudoscalar 
mesons one has to calculate eq. (122) and then solve eq. 
(119). The mass obtained by this procedure has not to 
be real. Indeed, when the mass of the meson is larger 



than that of its constituents the meson can decay into its 
constituents 

After carrying out the frequency sum v n ,q) can 

be brought as well in the form of eq. (115) with I\ and 
I 2 given by 



ii(ra) 
J 2 (rai,ra 2 ) 



d 3 p 1 
(2tt) 3 2Ep 



[l-f(E p -n)-f(E p + »)] 



d 3 p 



1 f(E p + fj) + f(E p - fj) - f(E p+q + fj) - f(E p+q - fj) 



(2tt) 3 2E p 2E p+q u - 

d 3 P l-f(E p -ri-f(E p + q + ») 



E, 



(2tt)* 



2E p 2E p+q 



1 



l€ 



(123) 



E n 



J P+q 



E n 



J p+q 



with E p = \/ml + p 2 and T^ +g = \/ m 2 + (P + ( l) 2 - I n 
the present approach we limit our mesons to the pseu- 
doscalar mesons. 

The model contains five parameters: the current mass 
of the light and strange quarks, the coupling constants 
Gd and Gs and the momentum cut-off A, are fixed by 
physical observables : the pion and kaon masses, the pion 
decay constant, the scalar quark condensate (qq) and the 
mass difference between 77 and 77' '. We will employ the 
parameters set : m q = 5.5 MeV, = 140.7 MeV , 
G s /A 2 = 1.835, G D /A 5 = 12.36, A = 602.3 MeV. The 
masses for up and strange quarks, as well as for ir and 
K for this parameter set [11, 32] are displayed in Fig. 
5. We see that for small \i and T the meson masses are 
smaller than the masses of their constituents. For large 
fx and T the opposite is true. If the mass of the con- 
stituents become smaller than the meson mass the me- 
son mass becomes complex and the meson become quasi 
particles. They exist is the plasma but with a life time 
which decreases with increasing \i and/or T (the width 
T = 2Gs^sTl p / s is displayed in yellow in Fig. 5). 



are given by 



(124) 



-iM t = ^ Ci ,c 3 ^c 2 ,cMps) Tu (Pi) 

+ 8 Cl ,c 3 Sc2,c4u(p 3 )(i'y 5 T)u(p 1 ) 
[iV[( Pl -p 3 )]v(p 4 )(i l5 T)v(p 2 ) 
-iM s = 5 CuC2 S C37C4 v(p 2 )Tu(p 1 ) 
\iV s s {p 1 ^p 2 )]v{p A )Tu{p 3 ) 
+ S Cuc2 S C3 ^v(p 2 )(ij 5 T)u(p 1 ) 
\ iV u(Pi + P2)]v(p4){i^T)u{pz), 



where pi(p 2 ) is the momentum of the incoming q(q) and 
^3(^4) that from the outgoing q(q). The are the color 
indices and T are the isospin projections on the mesons. 
V s and V p are the meson propagators of the form 



v s/i 



2G S 



1 - 2G s U p / s 



(125) 



C. Cross sections 

If created in heavy ion collisions the quark gluon 
plasma will expand rapidly. Therefore, not the static 
properties of the theory but the cross sections between 
constituents become dominant. In the NJL model these 
cross sections can be calculated via al/iV c expansion [35]. 
All the details can be found in [13, 35, 36]. Therefore we 
mention here only the essential facts. 



1. Elastic Collisions 

The Feynman diagrams for the qq — >• qq cross sections 
are displayed in Fig. 6. We see contributions from the 
s-channel and from the t-channel. The matrix elements 



with n p / 5 being the polarization tensor in the pseudo- 
scalar/scalar channel. This cross section is displayed in 
Fig. 7. We see that for most of the center of mass ener- 
gies this cross section is of the order of some mb. Close 
to the Mott transition the cross section increases dra- 
matically to more than hundred millibarn. The reason 
is the s-channel in which the incoming quarks become 
resonant with the intermediate meson [14]. This increase 
we observe in all channels, see Fig. 7. This means that 
at the end of the expansion of the plasma, shortly before 
the Mott temperature TMott is reached, the system comes 
almost certainly to a local equilibrium. Whether at tem- 
peratures much higher than TMott en local equilibrium 
can be established or maintained is in view of the size 
of the cross section not evident. The elastic qq and qq 
cross sections are of the order of a couple of mb. Because 
they do not have a s-channel they do not increase close 

to TMott 
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FIG. 6: The Feynman diagrams for elastic qq scattering. 




FIG. 7: Elastic cross section for the different channels as a 
function of y/s for a temperature close the Tjytott and at \i — 
[36]. 



2. Hadronization Cross Section 

The 1/N C expansion provides as well the hadronization 
cross sections in which a qq pair creates two pseudoscalar 
mesons. The Feynman diagrams are displayed in Fig. 8. 
For the details of the calculation we refer to [13, 37]. Fig. 
9 displays the cross section uu — » 7r + 7r - as a function 
of y/s for different temperatures. We observe that the 
cross section increases close to the kinematical threshold. 
Close the Mott transition the cross section can reach 100 
mb. Although the NJL model has no confinement this 
large cross section provokes that close to the cross over qq 





FIG. 8: The Feynman diagrams for hadronisation. 
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FIG. 9: Example of inelastic (right) [14] qq cross section as a 
function of Js. 



pairs create mesons very effectively (the cross section for 
the backward reaction is kinematically suppressed) and 
therefore most of the quarks are converted into mesons 
when the system reaches the Mott transition. 

We created tables of all elastic and inelastic cross sec- 
tions for our simulations. 
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IV. SIMULATION PROGRAM 

In this section we discuss how we use the results of 
the previous sections to formulate a molecular dynamics 
approach to describe the expanding q/q plasma. The 
expansion of such a plasma which is presumably created 
in the reaction between two heavy ions at ultrarelativisitc 
energies (at y/s well above 10 GeV ) is presently highly 
debated. 

Hydrodynamical calculations describe the properties 
of particles with low transverse momentum quite reason- 
able, at RHIC energies (yfs = 200 GeV) as well as at 
LHC energies (y/s = 2.46 TeV). There are, as discussed 
in the introduction, many flavours of hydrodynamical ap- 
proaches (ideal and viscous hydrodynamics, separation of 
core and corona, fluctuating or smoothed initial condi- 
tions, different equations of state close to hadronisation) 
which give very comparable results because different as- 
sumptions can have quite similar effects on observables. 
This makes it difficult to identify the physical processes 
from an agreement of the final result with the observ- 
ables. 

We describe the expansion in a relativistic quantum 
molecular dynamic approach, discussed in section 2, 
which is based on the NJL Lagrangian, discussed in 
section 3. We assume that the system remains suffi- 
ciently close to a local thermal equilibrium that we can 
parametrize the masses of the quarks and mesons by a 
local temperature and a local chemical potential. The 
quarks interact in two ways: Firstly, they change the 
mass of fellow quarks by their contribution to the chem- 
ical potential and to the local temperature, secondly, 
they interact with their fellow quarks by elastic and in- 
elastic scattering. This transport approach is called IN- 
TEGRAL (INTEractive Generalized Relativistic ALgo- 
rithms). 

The basic structure of a molecular dynamics program 
is described in the Fig. 10. We discuss in the following 
each of these steps. 



A. Initial conditions 

Principally any phase space distribution of partons can 
be taken as an initial condition for our calculations. In 
these first studies we assume a quite smooth initial dis- 
tribution which is determined as follows. In a first step 
we calculate the radius of the colliding (identical) nuclei 

by 



R 



ro 



A 1 ' 3 



(126) 



with ro = 1.25. For a finite impact parameter b we ap- 
proximate the overlap region by an ellipse with the axis 
yi E and y E : 



x# = R-b/2 



y E = ^(R-b/2)(R + b/2). 



(127) 



(initialization) 




(Analysis j 



FIG. 10: Standard algorithm for molecular dynamics calcula- 
tions. 



The extension in the third dimension is assumed to be 
proportional to the creation time of the QGP. We take 



z E = 2 tq = 2 fm. 



(128) 



Knowing the atomic number of the colliding nuclei, A, 
and the impact parameter, 6, we can construct the over- 
lap zone in coordinate space. 

In the present study we we assume that the system 
is close to local equilibrium. The mass of the partons 
and, consequently, their energy depends then on the lo- 
cal (T, /i), which depend on the initial To and /no for the 
centre of the collision. The initial local temperature de- 
pends on the position of the parton. At a point (x,y) 
with r = y/x 2 + y 2 the local temperature is a function of 
r/ro where the vector tq points in the direction of r and 
ro = \/ x i; + y%- The initial temperature is given by 



T(r) 



(129) 



l+exp(10 (r/ro-0.8))' 
The central initial temperature Tq can be parametrized 



as 



To(MeV) = 68 



log (Vi^(GeV) + 1) 
l + exp(1.5 (l-x E )) 



(130) 



Fig. 11 shows the initial temperature as a function of 
r/ro, Fig- 12 is a contour plot of the initial temperature. 
For RHIC energies this parametrization corresponds to 
the initial temperature of hydrodynamical calculations. 
Knowing the critical temperature T c = 165 MeV, this 
equation give us for RHIC (y/s^N = 200 AGeV) : Tq ~ 
2.2T C , and for LHC (^/siViv = 2760 AGeV) : T ~ 3.5T C . 
The calculations we present here are calculated with /i = 
H(r) = 0. 

Knowing (T(r), \i = 0) we can determine the mass and 
the initial density of quarks and antiquarks. 



N f 

= v =9 J 



fHP) 



d 3 p 



(2ir) 3 (hc) 3 



(131) 
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FIG. 11: Distribution of the temperature T as a function of 
the normalized radius r/ro. 




FIG. 12: Distribution of the temperature T (in MeV) in the 
transverse plane x,y. 



where f ± is the Fermi-Dirac distribution (eq. Ill), 

/ ± (p, m(T)) = 1 + exp ^Vp 2 + m2 =•= /V^) 7 

(132) 

g the degeneracy of the considered parton, TV the quark 
number and V the volume in the centre of mass system 
of the reaction (an ellipse with constant thickness) 



V = 7r x E y E z E . 



(133) 



Knowing density and volume of a slice with a given 
temperature we determine the local number of partons 
with help of eq. (129). The partons are then placed 
randomly in this slice. Each parton obtains its initial 
momentum p applying a Monte Carlo procedure which 
models the local Fermi-Dirac distribution with (T, /i). In 
the spirit of the core corona model close to the surface 
we assume thermalization only in longitudinal z direction 
and limit the transverse momentum in outward direction 



by limiting the azimuthal angle. This procedure en- 
sure that fast partons in corona are comovers and can 
hadronize easily. The spatial distribution of these par- 
tons is quite smooth. That is why we call this initial 
condition model the Hot Pancake Model (HPM). 



B. Transport model 

The partons in the expanding system are described 
by their positions and their momenta. The equations of 
motion of the particles are given by eq. (99) 



dr 
dr 



Ei 



N 

£ 

fe=i 



1 dV k (q' T ) 



(134) 



2 Eh 



In the NJL model the potential between the particles is 
a scalar interaction. This interaction acts like a mass 
which depends in our local equilibrium assumption on 
temperature and chemical potential of the environment. 
Therefore we can reformulate our energy constraints eq. 
(94) by 

Ki = - mf(T, M ) = 0. (135) 

This modifies the equation we have to solve : 



dpf 
dr 



N 

£ 

k=i 



m 



k dm k 



Ek dqi 



with 



dml _ dm% dT k 



dm* k dfjbk 



dqi 



dT k dqi dfJL k 



(136) 



(137) 



and with Tk{nk) being the local temperature (chemical 
potential) of the environment of the particle k. The de- 
pendence of the mass of the partons as well as that of 
it' s and K' s on the chemical potential and on the tem- 
perature is displayed in Fig. 13. Here we assume that 
the chemical potentials of up, down and strange quarks 
are identical. The masses show the expected behaviour 
of a cross over at high T and \i c± 0. 

At high (T, \i) we see the bare mass of the partons. 
When approaching low (T, \i) we observe a steep rise of 
the mass due to the scalar potential which becomes finite. 
At (T = 0, fi = 0) the light partons have their constituent 
mass of around 370 MeV. The 't Hooft term connects up 
and down quarks with strange quarks. Therefore the 
dependence of the strange quark mass on temperature 
and chemical potential becomes more complex. We see 
a first steep rise of the mass when the chemical potential 
arrives from above the transition temperature of the light 
quarks u : d and a second rise when the genuine transition 
of the s quark takes place. 

To solve the eqs. (99), the differential equations are 
converted in finite difference equations with a variable 
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FIG. 13: Dependence of the masses on the temperature and on the chemical potential of the environment in the NJL model. 
In the top row we display the masses of u and s quarks, in the bottom row that of ty and K mesons. 



time step. Its definition will be discussed in subsection 
IV E. For the solution we employ an adaptive method, 
depending on the time step size, with either a Euler al- 
gorithm or a Runge-Kutta algorithm of second (RK2) 
or of forth order (RK4). The cross sections and masses 
which have been calculated in [32] , [38] and [39] have been 
tabulated as a function of (T, /i, y^) and a linear inter- 
polation has been applied to accelerate the calculations. 



N 



d 3 p 



V * J (2nf(hcf 
2(f+(p,m u )-f-(p,m u )) ( 139 ) 

+ {f + (p,m s ) - f~(p,m a )) 



C. Thermodynamical medium 

In our local equilibrium approximation the effective 
mass m* of the partons depends on the temperature and 
chemical potential of the local environment. Therefore 
we have to construct these two quantities from the infor- 
mation on the system which is available, the 4-positions 
and 4-momenta of all particles. For this we define two 
densities, the fermionic density pp and the baryonic den- 
sity p B 



Pf{T,p) = — 



V 



d 3 p 



2(/ + (p,m n ) + r(p,m n )) ( 138 ) 



+ {f + (p,m 8 ) + / (p,m 8 )) 



with the degeneracy factor g = 2 x 3 = 6, and f ± defined 
in eq. (132). Neither pp nor pB are Lorentz invariants. 
In order to express Ti and pi as a function of the phase 
space coordinates (q^Pj) the following procedure is ap- 
plied : we introduce a Lorentz invariant Gaussian func- 
tion Rij(q f T ) inspired from the Wigner density eq. (8) 



exp 



<iTi 



L 2 



(140) 



to calculate the contribution of a neighbouring parton j 
to the density of the parton i. For the width we take 
L = 0.5 fm which is about the electromagnetic radius of 
known hadrons. This allows for a rewriting of the density 
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FIG. 14: Local temperature at r = as a function the dis- 
tance r and for a different number of fellow particles which 
have all a distance r to the considered particle. 
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FIG. 15: Collision between 2 particles in a medium. 



and therefore the mass has to be also updated in each 
time step for each routine (collision, decay and motion). 

How the distance between particles is related to the 
temperature can be demonstrated assuming that there 
are 1, 2 or 3 particles which have an identical distance 
r to the considered particle. This is shown in Fig. 14. 
The derivatives of the temperature with respect to the 
phase space variables, necessary to solve eqs. (134), are 
developed in detail in appendix 3. 



pFi — ^ Ri 



pBi = ^2 Ri i Si S n 0') > 

3& 



(141) 



with 



Sign(j) 



1 for fermions 
-1 for antifermions. 



(142) 



For /i = only one of these densities is necessary to 
determine the temperature. We use for this the Fermi 
density. Our approach corresponds to a Gaussian smear- 
ing of the density of a particle. These formula apply to 
free quarks and antiquarks. We also have to consider the 
partons which are bound in hadrons. For practical rea- 
sons, especially to avoid a sudden increase of the density 
and hence the temperature when mesons are produced, 
we consider mesons like one parton. 

Knowing ppi and pBi, eqs. (138) and (139) allows to 
determine Ti and pi . For p = and T » m the relation 
between Ti and pFi is analytical (see appendix 2) : 



PF= I 



9 



(143) 



where £ — 0.90154 being a normalization factor for the 
Fermi integral and the degeneracy factor becomes g = 
2x2x2x3x3 = 36. Then we find 



Ti = (he) 




(144) 



In the general case eqs. (138) and (139) have to be solved 
numerically. Ti and pi vary from time step to time step 



D. Cross sections and decays 

In addition to the potential interaction, which gener- 
ates the mass of the partons, the partons interact also by 
collisions. Collisions are characterized by cross sections. 
As in all transport theories these cross sections are con- 
verted into a geometrical concept which allows to decide 
which and when particles collide [16]. If two particles 
come closer than Ar = y 1 o ~/tt a collision between the 
particles takes place. In the program the collision is ex- 
ecuted at that time point at which the distance between 
the particles is minimal. 

In the present approach we have 4 types of processes : 

• q q -> q q, 

• q q^r qq, 

• q q — » M M (and backward process), 



M 



q q, 



where quarks are characterized by g, antiquarks by 
and mesons by M. These collisions increase the number 
of partons because for an expanding plasma qq^rMM 
, in which two partons are produced, is dominating over 
the backward reaction. Elastic collisions are primarily 
responsible for the thermalization of the plasma whereas 
the inelastic collisions are responsible for the hadroniza- 
tion. Both cross sections are small at temperatures well 
above the Mott temperature and therefore thermalization 
should happen only at the last stage of the expansion of 
the plasma shortly before the system hadronizes. 

Fig. 15 and 16 show a schema of a binary collision 
and a decay. The environment of both particles which 
enter a collision may be different and therefore we do not 
expect the same T and p for both particles. In order to 
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FIG. 16: Decay of a meson in a medium. 



determine the cross section which depends on T and \i 
we average over both particles. 

In the NJL approach quarks are not confined. Nothing 
prevents them from expanding into the vacuum. Nev- 
ertheless, applying our cross sections to the expanding 
system we find that at the end of the expansion almost 
all partons are bound in hadrons. The reason for this 
is the very large cross section for hadronization close to 
T c . Hence, when the system expands, close to T c hadron 
production becomes important. Hadrons formed slightly 
above T c live sufficiently long to survive until the system 
has passed T c and they become stable. 



E. Mean free path and time interval 

The geometrical interpretation of the cross section re- 
quires a careful study of the time step of the simulation. 
This can be seen by performing calculations in a box with 
periodic boundary conditions. As shown in Fig. 17 for 
the same initial condition (size of the box : a = 3 fm, 
filled with 30 free particles, for a duration of 10 fm/c) 
the total number of collisions depends on the chosen time 
step. On the top (bottom) we see the total number of col- 
lisions performed in the simulation program for the same 
initial condition as a function of the time step Ar and 
for a total cross section of 1 (5) mb. We miss collisions 
if the time step is above a critical value of Ar 

The reason for this observation is that if the mean free 
path is smaller than the time step, it is possible to have 
more than 1 collision per time step for the same particle, 
but numerically we only apply the first collision. There- 
fore the time step must be smaller than the mean free 
path. In our simulations the time between two collisions 
is given by the mean free path £ 

Ar = £ = (a p)' 1 , (145) 

which yields 

• Ar = 10 fm/c for a — 1 mb, 

• Ar = 2 fm/c for a = 5 mb. 

Fig. 17 show that in order to have the correct number 
of collisions the time step has to be much smaller than i. 
We need 

• Ar opt = 5.10 -2 fm for a = 1 mb , 
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FIG. 17: Total number of collisions for the same initial con- 
dition in a box simulation as a function of the time step Ar 
for a constant cross section of 1 mb (top) and 5 mb (bottom) . 



• Ar op t = 10 2 fm for a = 5 mb . 

To be on the safe side in our simulations we use 

Ar opt = (1000 (v reL ) (a) (p))' 1 (146) 

with the mean cross section (a) calculated from the col- 
lisions during the previous time step, the mean relative 
velocity (v re \.) and the mean scalar density (p) calculated 
for each time step. 

V. RESULTS 
A. Set up of the simulations 

The results we present here are obtained for simula- 
tions of Au-Au collisions at RHIC energies, ^snn = 
200 AGeV, or for Pb-Pb collisions at LHC energies, 
y's^jv = 2760 AGeV. We use the HPM (see subsection 
IV A) for the initial condition. 
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FIG. 18: Distribution of the mass (top left), of the momentum (top right), of the temperatures (bottom left) and of the scalar 
density (bottom right) for light and strange quarks as a function of the distance r from the centre and for central collisions. 



Fig. 18 displays results for central collisions. We show 
the mass, the momentum, the temperature and density 
of the light (ii, d) and heavy (s) quarks for RHIC and 
LHC initial conditions as a function of the position of the 
quarks measured with respect to the centre of the colli- 
sion. In the centre the quark mass is close to the bare 
mass. The more the surface is approached, where the 
density is smaller and the temperature lower, the mass 
increases and close to the surface we approach the con- 
stituent quark mass. The decreasing temperature is also 
responsible for the decrease of the average momentum. 
Finally, a smaller temperature yields also a decrease of 
the density when the surface is approached. 

In passing we notice that our smooth initial condition 
(Fig. 12) is not very realistic if one wants to compare 
quantitatively our results with experiments. Calculations 
in models like EPOS [40] show initial energy fluctuations, 
as seen in Fig. 19. Such fluctuations are visible in the 
final spectra of the mesons and can therefore not be ne- 



glected. We leave calculations with such realistic initial 
conditions to future investigations. 

Fig. 20 shows the number of initial partons and fi- 
nal particles (partons or hadrons) as a function of the 
impact parameter (the number of partons can increase 
due to the decay of mesons). The number of particles 
increases strongly with the centrality of the collisions 
and therefore also the computing time. The complex- 
ity of molecular dynamics calculations manifests itself in 
the fact that due to the two body potential and due to 
the collisions the computing time is quadratic in 0(N 2 ). 
Adaptive algorithms have been developed which reduces 
this dependence for large systems to 0(N log(iV)) [41] 
or even to O(N) (as the program of [42]). Unfortunately 
our approach does not fulfil the condition for such a sim- 
plification. Therefore central collisions are very time con- 
suming and only little statistics has been acquired up to 
now. 

One may use parallelization on supercomputers, or 
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more recently on graphic cards, to make the calculation 
feasible. In such an approach the calculation is decom- 
posed in small subtasks which can been executed inde- 
pendently. In our approach most of subtasks can be par- 
allelized and therefore the dependence of the calculation 
time on the particle number can be reduced by a factor 



R 



1 



(1-*) + 



(147) 



N 



where N is the number of processors and s the fraction 
of parallel tasks (in percent). In the ideal case of a code 
which is to 100% parallel and N is very large, the calcu- 
lation time is then O (TV 2 ) — » O(N). This parallelization 
requires the reprogrammation of parts of the program 
and is therefore a project for the future. 



FIG. 19: Energy density distribution (in GeV/fm 3 ) in the 
transverse plane x, y for EPOS [43] for a RHIC collision at 
6 = 6 fm. Coloured areas are QGP bubbles. 



4000 




FIG. 20: Number of initial and final particles (top) and com- 
putation time per simulation (bottom) as a function of the 
impact parameter b on a computer with monocore-CPU@3 
GHz. 



B. Check of the algorithm 

The most important check for the consistency of the 
derivation and its numerical realization is the energy con- 
servation. In molecular dynamics calculation it has to be 
strictly conserved. Fig. 21 displays the variation of the 
total energy of the system as a function of time for a 
simulation of a central RHIC collision. Such a simula- 
tion contains a couple of thousand partons. We see on 
the right plot that the energy varies by less than 0.2 %. 
The small variation of the total energy doesn't come from 
the solution of the differential equation (Euler or Runge- 
Kutta), but from the local density jump when a meson 
decay appears in a "low" density area. 



C. First results 

In this section we present some preliminary results 
which we have obtained for RHIC and LHC energies. 
They show that basic observables are well reproduced in 
our approach. In Fig. 22 we display the elliptic flow, 
i>2, as a function of the impact parameter. These re- 
sults are compared with the experimental data from the 
PHOBOS collaboration [44]. We see that the results of 
our approach agrees quantitatively quite well with the 
experimental finding. In this plot error bars are not the 
root mean square, but the variance of mean value af- 
ter N simulations. This explains the huge error bars for 
small impact parameters because of the large number of 
particles and the small number of simulations N. 

Fig. 23 displays how the elliptic flow develops as a 
function of time. On the top we display our results and 
on the bottom that one of PHSD calculations [45]. By 
definition there is initially no elliptic flow (assuming ther- 
mal equilibrium). In both calculations the flow develops 
very similarly and both calculations agree also on the fi- 
nal value. In our case, despite of the small cross section 
of about 1-2 mb, we observe initially many collisions due 
to the high density. They thermalize the plasma rapidly 
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FIG. 21: Variation of the total energy for a central RHIC 
collision as a function of time. 
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FIG. 23: Time evolution of the elliptic flow V2 for b — 6 fm 
(top), and the comparison with PHSD calculations for similar 
initial conditions [45] (bottom). 



FIG. 22: V2 as compared to experimental data from the PHO- 
BOS experiment [44] as a function of the impact parameter 
b. 



and lead to an elliptic flow in less than 1 fm/c. The flow 
is lowered later by the change of the NJL masses. 

Experimentally it has been found that the transverse 
momentum spectra of 7r and K have a different shape 
[46]. This can be seen in Fig. 24, bottom, where we 
compare the experimental data with results from hydro- 
dynamical calculations [47]. On the top we display our 
results. We observe the same difference of the slopes 
as seen in experiment which is usually attributed to the 
hydrodynamical evolution of the system. The lines are 
there to guide the eyes. 

Fig. 25 displays a contour plot of the number of col- 
lisions as a function of the distance to the centre of the 
initial ellipse r and of the time t. On the top we display 
inelastic collisions on the bottom elastic collisions. 

Initially we have a very high density zone where elastic 
and inelastic collisions take place frequently because the 
mean free path is small despite of the small cross section. 



When the system expands the density becomes lower but 
the cross section does not increase. Therefore we observe 
less collisions. When the system approaches the critical 
temperature the cross sections becomes very large, this 
compensates by far the decrease of the density and there 
the collision rate becomes large again for elastic as well 
as for inelastic collisions. Here the hadrons are created 
which finally survive. 

For the LHC initial condition, Fig. 26, we see the 
same phenomenon but a longer lifetime of the plasma. 
In contradistinction to simulation for RHIC energies the 
corona partons do not hadronize early but the stream of 
partons from the interior heats up the surface. So the 
system expands in the quark phase and hadronization 
takes place only much later over a large space time area. 

Fig. 27 shows the distribution of y/s and of T at which 
the final hadrons are produced. We see a broad distri- 
bution around the critical temperature and not a single 
freeze out temperature as assumed in the Cooper-Frye 
formula [48] which is used to created hadrons in hydro- 
dynamical calculations [47]. The temperature at the K 
production points is slightly lower than that for the tt 
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FIG. 24: The diY/27rp T dp T spectrum for b = 4 fm (top), 
and the results of a hydrodynamical calculation and of the 
experiment for similar conditions (centrality 0-5%) (bottom) 
[47]. 



production points as expected by the NJL cross sections. 

Fig. 28 compares the hadronization in our approach 
(top) with the results of PHSD calculations (bottom) 
[45]. We observe in both calculations a hadronization 
time of around 5 fm/c for the RHIC conditions. The dif- 
ference between the particle number of PHSD and in our 
model comes from the fact that we have a different ini- 
tial conditions. Of course the hadronisation is longer for 
the LHC initial condition (top right) due to the higher 
density of partons. 



VI. SUMMARY 

We have presented in this paper a relativistic molecu- 
lar dynamics approach. We have shown that for a spe- 
cific choice of the constraints it is possible to find back 
the classical relativistic equations of motion. These con- 
straints give us physical trajectories with causal motion 
and the conservation of the energy of a strongly interact- 
ing system. Using the Nambu-Jona-Lasinio Lagrangian 
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FIG. 27: y/s distribution for inelastic collisions (top), and the 
distribution of the temperature at the production points, T, 
for pions and kaons (bottom). 



to describe the potential interactions and the scattering 
among the partons we find that is it is possible to model 
the expansion of a quark anti-quark plasma. Close to the 
cross over the elastic as well as the hadronization cross 
sections explode. The large hadronization cross section 
is the reason that the large majority of the quarks form 
mesons which can finally be observed. 



We presented first result which show that a approach 
which does not enforce local equilibrium, as hydrodynam- 
ics, and in which the transition to the hadronic world is 
not sudden, as in the Cooper Frye approach, used fre- 
quently in hydrodynamical calculations, gives qualitative 
agreement with some key observables. Further studies for 
a in depth comparison with the existing models will be 
the subject of a future publication. 



The approach is in spirit close to the PHSD approach 
but differs completely as far as the temperature and den- 
sity dependence of the mass of the quarks is concerned. 
Therefore it will be fruitful to compare the observables 
which are obtained in both approaches for the same ini- 
tial condition. 
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Appendix 



1. Relativist ic calculations 

a. Derivatives of transverse distances 

The calculation of the derivatives for transverse dis- 
tances can be done rigorously 



dq 2 



dqkv 



= 2q 2 



= 2q 2 



= 2qi 



dqi 



dqkv 



dqj 



'Ikv 



dqkv 



(ft/ ~ (Qij^jKj) 



d (qijo 



dqkv 



dqkv 



(148) 



dq' T 



dpkv 



dq' T 



dpkv 
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M dp^ 
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d(q i3 aUn 
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= -2q: 



/ v (qijaU a ) lab 



Tij ' 



= 0. 



b. Matrix of constraints 



(151) 
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dp, 
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(149) 



We present the full calculation of the matrix of con- 
straints and the full expression of the equations of mo- 
tion for the case in which the KT [26] condition is not 
fulfilled. We start with the calculation of the derivative 
of the first constraints 



dK % 



dKi 
dPk 



dPk 



Then for time constraint 



(152) 



and the same kind of derivatives can be found for q' T : 
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and 



dXi V^/r 

d X i 1 



ik 
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dpi NVP 2 z 



dxN = Up_ 
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dXN 1 ST v Q 
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We notice that except dxi/dq%, the derivatives of x don't 
depend on k. For the full matrix of constraints we find 



{K uXj } = J2 

k L 



jin 



N 



v ^ 1 



E ^j/p 



=0 



^ dpi J N 



! *E(ii,-4)l. 
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(y dW] Up 



N 



N 



Using eq. (82), we can write : 
{Xi,Xj} = E 

k 

(E(^-^)^r^f E^ 



(155) 



(156) 



(157) 



{Xi,Xn} = 


E 


( 





(158) 



1 V ,. 
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(159) 



We can summarize these results presenting the com- 
plete matrix of constraint : 



{Ki, Kj} {Ki, Xj } 
{Xi,Kj} {Xi,Xj} 



with 



A- 1 ={K i ,K j }^Q 
B-/ ={Xi,Xj} = . 



(160) 



(161) 



S^ = { Xi , Kj }^0 
That gives us the following relativistic factor : 

A fc = C 2N k K k < 27V. (162) 
The final expression for the equations of motion is : 
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We can write these equations in a simplified form us- 
ing q' T and the global reference frame (f7 M = (1,0,0,0)). 
Please notice that we only use the 3- vector part of these 
equations : 



M 

dr 



2N-1 

E^J E A * 
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(164) 



E 

fc=l 



1 dV k 
2E k dq iu 



The second term of dgf/dr is embarrassing. For N = 2 
particles this term is vanishing. Unfortunatelly for a large 
number of particles (2 < N) it does not disappear. To 
avoid this the KT condition must be fulfilled. 
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2. Thermodynamical densities 



3. Derivation of potential 



Eq. (138) can be calculated analytically for \i — > and 
m <C T (including a factor 2 in g) : 
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Then we find : 
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For the baryonic density, eq. (139), we have 

47T 



(27r) 3 (ftc> 



/<OC 

,5/ (/ + -r)p 2 d P 

./o 



4tt tT 3 ( fl 



47T 7T 2 



/i\3 1 



(2^) 3 (ftc) 3 * 3 
^6tt 2 VfiJ ' 



Low Order 



and finally 



\Li (He) 



r 2\ 1/3 







(165) 



(166) 



(167) 



(168) 



The degeneracy factor for the spin, the parity, the color 
and the flavour is # = 2x2x3x3 = 36. 



The forces coming from the derivative of the potential 
are : 
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The local temperature is defined as : 
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with Rij = exp(q f T ' 2 ./L 2 ) and 
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which can be rewritten 
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n=(ac)f— J ~190MeV. (171) 

L is a weighting factor. Its derivatives are : 
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